/**
 * @file	IIRFilter.h
 * @author  Thomas Fu <thomas.ks.fu@gmail.com>
 * @version 1.0
 *
 * @section LICENSE
 *
 * This program is propety of The SOSS and may not be used by any 
 * individual for any purpose at any time without express permission from 
 * Mike Casey, Thomas Fu, Chris Hairfield, Kyle Lamson, Ben Treweek, or Cliff 
 * Warren.
 *
 * @brief
 * This class is a representation of a digital, infinite impulse response (IIR)
 * filter. 
 *
 * @class IIRFilter
 *
 * @todo
 * Fix the representation of polynomial coefficients across functions to be 
 * consistent (either all double*'s or all vectors)
 */

#ifndef _IIRFilter_H
#define _IIRFilter_H

#include <cstdlib>
#include <vector>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>

using namespace Eigen;

using namespace std;

/** Enumerated type used to specify the type of filter to be constructed */
enum FilterType {LOWPASS, HIGHPASS, BANDPASS, BANDSTOP};

class IIRFilter
{
	public:
		
		/**
		 * Constructor for the IIR Filter class.
		 * @param num A vector of doubles containing the coefficients of the 
		 * polynomial numerator of the z-domain transfer function of the IIR filter.
		 * 
		 * @param den A vector of doubles containing the coefficients of the 
		 * polynomial denominator of the z-domain transfer function of the IIR 
		 * filter. 
		 */
		IIRFilter(vector<double> num, vector<double> den);
		
		/** Destructor for the IIR filter */
		~IIRFilter();
	
		/**
		 * Computes filtered value of the given data point based on the values of 
		 * the filter's transfer function as well as a time history maintained
		 * internally by the filter. This means that tickFilter must be called on
		 * the samples in temporal order.
		 *
		 * @param val The next value in the signal to be filtered
		 *
		 * @return The filtered value corresponding to the unfiltered parameter 
		 * passed into the function
		 */
		double tickFilter(double val);
		
		/**
		 * Computes the coefficients of the numerator and denominator polynomials 
		 * for the discrete-time (z-domain) transfer function for a Butterworth 
		 * filter of arbitrary order.
		 *
		 * @param order The number of poles to be used in the filter (the greater 
		 * the number of poles, the sharper the filter response).
		 *
		 * @param w_c The cutoff frequency of the Butterworth filter
		 *
		 * @param w_s The sampling frequency of the signal being filtered by the 
		 * Butterworth filter
		 *
		 * @param type A value specifying the type of filter to be implemented
		 *
		 * @param num An array of doubles that will contain the coefficients of the 
		 * polynomial numerator of the transfer function after this function finshes
		 * executing. Space for this array must be allocated before this function is
		 * called.
		 *
		 * @param den An array of doubles the will contain the coefficients of the 
		 * polynomial denominator of the transfer function after this function 
		 * finishes executing. Space for this array must be allocated before this 
		 * function is called.
		 *
		 * @todo 
		 * Add funcationality for Bandpass and Bandstop modes
		 */
		static void ButterworthFilter(int order, double w_c, double w_s,
													 FilterType type, double *num, double *den);
	
		/**
		 * Computes the coefficients of the numerator and denominator polynomials 
		 * for the discrete-time (z-domain) transfer function for a first order 
		 * shelving filter.
		 *
		 * @param b Vector of doubles representing the coefficients of the numerator
		 * polynomial for the discrete-time (z-domain) transfer function
		 *
		 * @param a Vector of doubles representing the coefficients of the 
		 * denominator polynomial for the discrete-time (z-domain) transfer function
		 *
		 * @param tone A double between 0 and 100 that changes the output tone of 
		 * the shelving filter (0 corresponds to a more muddy, bassy sound whereas
		 * 100 corresponds to a more trebly sound).
		 *
		 * @param w_c The cutoff frequency of the shelving filter
		 *
		 * @param w_s The sampling frequency of the shelving filter
		 *
		 * @todo
		 * Add support for higher order shelving filters
		 */
		static void shelvingFilter(vector<double> &b, vector<double> &a, double tone, 
															 double w_c, double w_s);
		
	private:
		
		// Maintains the order of the IIR filter (one less than the length of the 
		// degree of the numerator and denominator polynomials)
		int order;
		
		// Maintains the index in xHistory and yHistory corresponding to the oldest
		// value
		int oldestBufferEntry;
		
		// Stores the coefficients of the numerator of the z-domain transfer
		// function
		vector<double> num;
	
		// Stores the coefficients of the denominator of the z-domain transfer
		// function
		vector<double> den;
	
		// Array of doubles used to store the time history of the input. For 
		// efficiency, this is maintained as a circular buffer, with the oldest
		// value always overwritten by the newest
		double *xHistory;
	
		// Array of doubles used to store the time history of the output. For 
		// efficiency, this is maintained as a circular buffer, with the oldest
		// value always overwritten by the newest
		double *yHistory;
	
		/**
		 * Generates all combinations of poles that occur in the generation of a 
		 * polynomial from a list of zeros. Effectively generates a power set of the
		 * set {1, 2, ..., n} where n is the number of zeros of the polynomial. This 
		 * can be used to generate a polynomial by recognizing that, for each term, 
		 * the power of the variable plus the number of zeros involved in the 
		 * coefficient must be equal to n, the total number of zeros. For example,
		 * for the polynomial (x + 1) * (x + 2) * (x + 3), one term will be 2 * x^2
		 * and another term will be 2 * 3 * x.
		 *
		 *
		 * @param length The number of zeros of the polynomial (including algebraic
		 * multiplicity).
		 *
		 * @return Returns a vector of (mulitplicative) combinations of zeros. Each
		 * combination is represented as a vector that lists the indices of each 
		 * zero to be multiplied for a given term.
		 */
		static vector< vector<int> > combinations(int numZeros);
	
		/**
		 * Converts a list of complex zero locations for a polynomial into an 
		 * array of doubles representing the coefficients of each term of the 
		 * polynomial when multiplied out. That is, it converts a polynomial from
		 * the form (x + 1)(x + 2) into the form x^2 + 3x + 2. 
		 *
		 * @param zeros A vector of complex<double>'s indicating the locations of 
		 * the zeros of the polynomial. Because this function is used to generate
		 * filters for physical systems, it is assumed that complex zeros always 
		 * come in complex conjugate pairs.
		 *
		 * @return An array of doubles indicating the coefficients of the polynomial
		 * corresponding to the zeros specified. The indexing convention is such 
		 * that the 0th index of the array returned corresponds to the coefficient
		 * of the highest power of the variable.
		 */
		static double* polynomial(VectorXcd zeros);
	
		/**
		 * Converts a list of complex zero locations for a polynomial into an 
		 * array of doubles representing the coefficients of each term of the 
		 * polynomial when multiplied out. That is, it converts a polynomial from
		 * the form (x + 1)(x + 2) into the form x^2 + 3x + 2. 
		 *
		 * @param zeros An array of complex<double>'s indicating the locations of 
		 * the zeros of the polynomial. Because this function is used to generate
		 * filters for physical systems, it is assumed that complex zeros always 
		 * come in complex conjugate pairs.
		 *
		 * @param numZeros An integer representing the number of zeros in the 
		 * polynomial
		 *
		 * @return An array of doubles indicating the coefficients of the polynomial
		 * corresponding to the zeros specified. The indexing convention is such 
		 * that the 0th index of the array returned corresponds to the coefficient
		 * of the highest power of the variable.
		 */
		static double* polynomial(complex<double> *zeros, int numZeros);
};

#endif